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Global Radiation-Magnetohydrodynamic Simulations of Black Hole 
Accretion Flow and Outflow: Unifled Model of Three States 
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Abstract 

Black- hole accretion systems are known to possess several distinct modes (or spectral states), such 
as low/hard state, high/soft state, and so on. Since the dynamics of the corresponding flows is distinct, 
theoretical models were separately discussed for each state. We here propose a unified model based on 
our new, global, two-dimensional radiation-magnetohydrodynamic simulations. By controlling a density 
normalization we could for the first time reproduce three distinct modes of accretion flow and outflow 
with one numerical code. When the density is large (model A), a geometrically thick, very luminous disk 
forms, in which photon trapping takes place. When the density is moderate (model B), the accreting gas 
can effectively cool by emitting radiation, thus generating a thin disk, i.e., the soft-state disk. When the 
density is too low for radiative cooling to be important (model C), a disk becomes hot, thick, and faint; 
i.e., the hard-state disk. The magnetic energy is amplified within the disk up to about twice, 30%, and 
20% of the gas energy in models A, B, and C, respectively. Notably, the disk outflows with helical magnetic 
fields, which are driven either by radiation pressure force or magnetic pressure force, are ubiquitous in any 
accretion modes. Finally, our simulations are consistent with the phenomenological a- viscosity prescription, 
that is, the disk viscosity is proportional to the pressure. 
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drodynamics: MHD — radiative transfer 



1. Introduction 

The extensive study of disk accretion flows started in 
the 1960's. The standard disk model, and then the slim 
disk model and the radiatively inefficient accretion flow 
(RIAF) model were proposed for explaining a variety 
of accretion modes (Shakura & Sunyaev 1973; Ichimaru 
1977; Rees et al. 1982; Abramowicz et al. 1988; Narayan 
& Yi 1994). These models are successful, but have some 
limitations. For example, the disk viscosity, the most im- 
portant key ingredient for the accretion disk theory, is pre- 
scribed by a phenomenological a- viscosity model, whereby 
the viscous torque is proportional to the pressure (Kato 
et al. 2008), although its physical basis is not clear. They 
are (radially) one-dimensional models so that they cannot 
describe multi-dimensional motion, such as outflow and 
internal circulation. Complex coupling between radiation, 
magnetic fields, and matters is not accurately solved, ei- 
ther. 

Since the disk viscosity is likely to be of mag- 
netic origin (Balbus & Hawley 1991), multi-dimensional 
global magneto-hydrodynamics (MHD) simulations are 
being rather extensively performed recently as a model 
for the disks with low luminosities (Matsumoto 1999; 
Machida et al. 2000; Hawley & Krolik 2001; Koide 
et al. 2001; De Villiers et al. 2003; Hawley & Krolik 



2006) . Such non-radiative MHD simulations cannot 
explain higher luminosity states, however, since strong 
matter-radiation coupling is expected. As an indepen- 
dent approach several groups performed two-dimensional 
radiation-hydrodynamic (RHD) simulations of very lumi- 
nous flow since the 1980's (Eggum et al. 1988; Okuda & 
Fujita 2000; Ohsuga et al. 2005; Ohsuga 2006). Those sim- 
ulations were, however, non-MHD simulations and so they 
were obliged to rely on the phenomenological a- viscosity 
model. Multi-dimensional radiation-MHD (RMHD) simu- 
lations are unavoidable. Such simulations were attempted 
in the past (e.g.. Turner et al. 2003; Hirose et al. 2006), but 
these are restricted to local simulations performed under 
the shearing-box approximations and, hence, global cou- 
pling of magnetic fields was artificially quenched there. 

We, here, report for the first time the results of global 
two-dimensional RMHD simulations with a motivation to 
establish a unified view of the accretion flow and outflow 
around the black holes. 

2. Numerical Method 

Our method of calculations is extension of that of MHD 
simulations (e.g., Kato et al. 2004). We use cylindri- 
cal coordinates (r, (/:?, z), where r is the radial distance, 
(f is the azimuthal angle, and z is the vertical distance. 
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We assume that the flow is non-self-gravitating, reflection 
symmetric relative to equatorial plane, and axisymmetric 
with respect to the rotation axis. General relativistic ef- 
fects are incorporated by the pseudo-Newtonian potential 
(Paczynsky & Wiita 1980). For the opacity, we consider 
the Thomson scattering, free-free absorption, and bound- 
free absorption (Rybicki & Lightman 1979; Hayashi et al. 
1962). The energy equations of gas and radiation are given 
by 



dt 



V • (^gas^^) 



An 



= -PgasV • V - AtTkB + CKErad + -^^J i (1) 



and 

dErad 
dt 



V • (^rad^) 

= -V • Frad - Vi; : Prad + ^TTf^B - Cf^E.^d. (2) 



where £^gas is the internal energy density of the gas, v is 
the velocity, Pgas is the gas pressure (= 2£^gas/3), B is the 
blackbody intensity, J is the electric current, £^rad is the 
radiation energy density, F^ad is the radiative flux, Prad is 
the radiation pressure tensor, is the absorption opacity. 
We adopt the anomalous resistivity, r^, which is the same 
as that used in Kato et al. (2004); 
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where Vd = J/p is the electron drift velocity, 'Ucrit = 0.01c 
is the critical velocity, and r^max = lO~^cRs is the maxi- 
mum resistivity with Rs being the Schwarzschild radius. 
This form of the anomalous resistivity was proposed by 
Yokoyama & Shibata (1994) to account for the occur- 
rence of fast reconnections in solar flares. The MHD 
related terms are solved by the modified Lax-Wendroff 
scheme (Rubin & Burstein 1967). We employ the flux- 
limited diffusion (FLD) approximation to solve the ra- 
diation energy equation (Levermore & Pomraning 1981). 
The radiation energy transport via the radiative flux is 
solved based on the implicit method, where we separately 
treat radiative fluxes in the radial and vertical directions 
with using Thomas method for a matrix inversion. The 
gas-radiation interaction is also solved with the implicit 
method, which is basically the same as that described by 
Turner & Stone (2001). An advection term in the en- 
ergy equation of the radiation is solved with the explicit 
method, in which an integral formulation is used to gen- 
erate a conservative differencing scheme. We performed 
the test of two-dimensional radiation propagation (Turner 
& Stone 2001), finding that the energy loss is less than 
0.005% in 600 steps, from which we estimate the error in 
the energy conservation not to exceed 0.08% during the 
photon traveling timescale in our simulations. 

The grid extends from 3Rs to lOSi^s- in the radial direc- 
tion and from to 91. 6Rs in the vertical direction. The 
grid spacing is uniform, 0.2Rs, in both directions. We 



adopt free boundary conditions for the matter and mag- 
netic fields; i.e., the matter can freely go out but not to 
come in and the magnetic fileds do not change across the 
boundary. We assume that radiation goes out with the 
radiation flux of c^^rad, except at r = 3Rs and z > 3Rs^ 
through which no radiation goes out. We assume a black 
hole mass to be IOMq. In future we will perform simu- 
lations, in which we will put the boundary condition at 
around Rs and to solve the both side of the midplane. 

We start calculations with a rotating torus, in which the 
magnetic fileds are purely poloidal (plasma-/? = 100) and 
closed loops in the torus, being located at around AORs 
embedded in non-rotating isothermal corona. Our initial 
conditions are the same as those of model B in Kato et al. 
(2004), except that the density of the corona is 0.05 times 
as large as that one. We evolve the initial torus by solving 
non-radiative MHD equations for 1 sec. We then assign 
the density normalization (po), density at the center of the 
initial torus, and turn on the radiation terms. We calcu- 
late three models in total, by setting po = lgcm~^ (model 
A), 10~^gcm~^ (model B), and 10~^gcm~^ (model C). 
Since radiation loss rate depends on the density, we can 
reproduce the three distinct regimes of accretion flow. 

3. RESULTS 

3.1. Accretion Flows 

Figure 1 clearly visualizes that the flow patterns differ 
significantly among three models. The typical mass ac- 
cretion rate (Mace), the luminosity (L), and the density 
(p) and the temperature (T), as well as other important 
quantities, are summarized in Table 1. 

In model A with a relatively large density normaliza- 
tion, the mass accretion rate exceeds the Eddington rate, 
Le/c^, with Le being the Eddington luminosity. The disk 
is optically and geometrically thick. Photons are not easy 
to go out from the surface due to a large optical depth 
so that radiative cooling is restricted. We confirm the 
photon- trapping effects. The disk is supported by radi- 
ation pressure. Circular motion appears in the disk re- 
gion. Since L>Le, this model corresponds to the two- 
dimensional version of the slim disk model. The calcu- 
lated temperature and density are also consistent. 

In model B with a moderate density normalization, a 
geometrically thin disk forms because of efficient radiative 
cooling. The disk is optically thick and supported mainly 
by radiation pressure, which is slightly greater than the 
gas pressures. Such properties, as well as the temperature 
and density, agree with those of the standard disk model. 
It might be noted well that the flow in model B has not 
reached a quasi-steady state, since the viscous timescale 
is about 45 sec at r = lORs, whereas the elapsed time is 
8 sec. In a forthcoming paper we will present finer mesh 
calculations with adequate grid spacing, by which we will 
be able to investigate the detailed, internal structure of 
the thin disk. 

In model C with a small density normalization, the den- 
sity is too low for radiative cooling to be important. The 
disk is filled with hot rarefied plasmas and is geometrically 
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Model A Model B Model C 




Fig. 1. Perspective view of inflow and outflow patterns near the black hole for models A, B, and C, from left to right, respectively. 
Upper panels: Normalized density distributions (color) are overlaid with isosurfaces, at which the outward velocity equals to the 
escape velocity. Here, the elapse times are 5.25 sec, 7.75 sec, and 5.0 sec for models A, B and C, respectively. Lower panels: The 
distributions of the normalized radiation energy density (color) is overlaid with the magnetic fleld lines. In models A and B, the 
green surfaces of the disks correspond to the photosphere, where the optical thickness measured from the upper (z = 91.6Rs) or 
lower (z = —91.6Rs) boundaries is unity. Note that the disk in model C is optically thin. 



thick but optically thin. We find significant circular mo- 
tion inside the disk. This model corresponds to the RIAF 
model. 

3.2. Outflows 

As shown in Figure 1, the disk outflows with helical 
magnetic fields are ubiquitous around the black holes. In 
models A and C the magnetic field lines stretch out verti- 
cally in the vicinity of the rotation axis. While, around the 
equatorial plane, the toroidal component of the magnetic 
fields is dominant over other components in all models, 
which are reminiscent of magnetic-tower jets (Lynden-Bell 
1996; Kato et al. 2004). 

In model A, the strong radiation pressure force is re- 
sponsible for driving the quasi-steady outflows above and 
below the disk, whose velocity amounts to ~ 0.25c. We 



find that the radiation energy density (£^rad) is very large 
in the disk region, and the steep profile of E^ad enhances 
the radiation force (radiative flux) (Ohsuga et al. 2005). 

Remarkably, our simulation of model B shows the occur- 
rence of magnetically powered disk wind, on the contrary 
to the usual belief regarding the standard disk model. 
Note, however, that the disk wind is not so strong in model 
B; Lkin <^ and Lkin/^ is the smallest among all models, 
where Lkin is kinetic luminosity. It must be stressed that 
we have solved the entire inflow-outflow structure simul- 
taneously unlike the previous simulations, in which the 
inflow (disk) structure was not solved but treated as the 
boundary condition (Proga & Kallman 2004). 

Our RMHD simulations reveal Lkin/^ > 1 in model C 
in contrast with models A and B, implying that the disks 
with Mace ^ Le/c^ lose the energy via the jets rather 
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Table 1. Three calculated models and basic properties. 



Model 




A 


B 


C 


Corresponding States 




Slim disk 


Standard disk 


RIAF 


normalized mass accretion rate* 


Macc/(i^E/c^) 


120 


5.3 X 10-^ 


5.8 X 10-^ 


Eddington ratio* 




1.0 


7.3 X 10-"^ 


7.4 X 10-^0 


ratio of outflow rate to accretion rate* 


Mout/Macc 


0.045 


0.013 


0.1 


ratio of kinetic luminosity to luminosity* 




0.16 


0.0013 


120 


density''' 


log/) [g cm~^] 


-1--2 


--5 


- -10 


gas temperature''' 


\ogT[K] 


-8 


-6 


10-11 


gas energy density-1- 


logEgas [ergcm"^] 


15.0 


9.7 


9.0 


magnetic energy density-'- 


log£;mag [ergcm"^] 


15.3 


9.2 


8.3 


radiation energy density-'- 


log^^rad [ergcm"^] 


17.2 


10.4 


3.2 


plasma- /3-'- 


Pgas / Pmag 


0.33 


2.1 


3.3 



* The mass accretion rate (Mace) is evaluated the mass passing through the inner boundary near the black hole (r = ^Rs and 
z < SRs) per unit time, and outflow rates (Mout) indicate the mass passing through the upper boundary (z = 91.6Rs) with higher 
velocities than the escape velocity (high-velocity outflow) per unit time. The luminosity (L) is calculated based on the radiative 
flux at the upper boundary, and the kinetic luminosity (L^in) is calculated from the amount of the kinetic energy carried by the 
high- velocity outflow. These values are time-averaged over t = 5.5 — 6.0 sec (models A and C) or t = 7.5 — 8.0 sec (model B). 
^ Time-averaged over t = 5.0 — 6.0 sec (models A and C) or t = 7.0 — 8.0 sec (model B) at the regions r = 10 — 20Rs and z ~ 0. 
■1- Density-weighted spatial averages within the regions of r = 10 — 20Rs and z = — lORs and are time-averaged over t = 5.0 — 6.0 
sec (models A and C) or t = 7.0 — 8.0 sec (model B). Here, pgas and pmag are gas and magnetic pressures. 



than via radiation. The outflow rate is 10 % of the mass 
accretion rate, and the ratio is largest in three models. 
The photons freely escape from the disk, producing quasi- 
spherical distribution of £^rad, whereas the radiation en- 
ergy is enhanced inside the disk in models A and B. The 
radiation force is negligible because of small radiative flux. 

The FLD approximation is good for Models A and C, 
since the whole region (except for the very vicinity of the 
inner boundary) is optically thick for Thomson scattering 
in the former and radiative cooling is never important in 
the latter. It is known to be problematic to determine the 
direction of the radiation flux in regions where the optical 
depth is around unity (e.g., around the disk surface in 
model B). We, however, wish to stress that the outflow is 
accelerated by the magnetic pressure, and not by radiation 
force. 

3.3. Amplification of magnetic fields and viscosity 

We find that magnetic energy (£^mag) is amplified to be 
30% and 20 % of the gas energy (£^gas) in models B and C, 
respectively (see Table 1). In model A, surprisingly, £^mag 
does exceed the gas energy, 



mag 



^2£^gas. Its implication 
is enormous: the viscosity had better been scaled in terms 
of the total (or radiation) pressure and not of the gas pres- 
sure in radiation pressure-supported disks (cf., Sakimoto 
& Coroniti 1981). 

It is one of the most significant issues in astrophysics 
how to prescribe the disk viscosity. In Figure 2 we show 
how the magnetic torque, {—BrB^/Air), behaves as a func- 
tion of the pressure in the region close to the black hole 
(at r = 5Rs)- We can see that the torque is roughly 
proportional to the total pressure (with some scatters) 
in all the cases. If we define the viscosity parameter as 
a = {—BrB^/4:7r)/{ptot) where total pressure (ptot) is the 
sum of gas pressure and radiation energy density divided 
by 3, we estimate a ~ 0.004, 0.006, and 0.002 for models 
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Fig. 2. The density- weighted spatial averages of the mag- 
netic torque near the black hole, r = 5Rs, as functions of the 
total pressure. Each point represents the time- aver aged value 
for every 0.25 sec. At the initial points, (\og[— {Br B^/Att)], 
log(j9tot>) ~ (15.3, 12.7), (10.9, 8.7), and (8.6, 5.8) in CGS 
unit for models A, B, and C, respectively. 

A, B, and C, respectively. 

Cautions should be taken to the point that the proper- 
ties of magnetic fields in the two-dimensional simulations 
could be affected by artificial occurrence of channel modes 
of flow and the anti-dynamo theorem and thus deviate 
from those of the three-dimensional simulations. Hence, 
the three-dimensional study should be explored in future 
work. 
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4. Observational Implications 

Outflows and jets are ubiquitous in any regimes of 
black-hole disk accretion flow. They seem to manifest 
themselves by warm absorptions features of active galac- 
tic nuclei and by blue-shifted absorption lines in the black 
hole binaries (BHBs) (Blustin et al. 2005; Cappi 2006; 
Kubota et al. 2007). The high- velocity outflows with ve- 
locity of ~ 0.25c in model A will explain the X-ray observa- 
tions of bright quasars exhibiting blue-shifted absorption 
lines, which are interpreted as the absorption by the out- 
flow material moving with velocity on the order of 0.1c 
(Pounds et al. 2003; Reeves et al. 2003). Although no 
strong outflows were expected in the framework of the 
standard disk model, Miller et al. (2006) concluded by X- 
ray observations of BHB, GRO J1655-40, that the X-ray 
absorbing wind is ejected from the geometrically thin disk 
viewed at an inclination angle of ~ 70°. The density and 
the velocity of the wind were reported to be ~ 10~^gcm~^ 
and 0.001c — 0.05c. Such features are roughly consistent 
with our simulation (model B), in which the gas with 
p = several x 10~^^gcm~^ is blown away towards the di- 
agonal direction (the polar angle of ~ 45°) at the speed 
of ~ 0.01c. Our simulations show that the magnetic field 
lines are along the jet axis (models A and C). Such struc- 
ture is revealed by the recent radio observations of polar- 
ized emission in BL Lac object, Markarian 501 (Giroletti 
et al. 2008). The detailed comparison with the observa- 
tions is left as future work. 

To summarize, our global RMHD simulations can open 
a new era of accretion disk research and provide a unified 
view of accretion flows in various contexts. 
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